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In the past decade, synchronization on complex networks has attracted increasing attentions from 
various research disciphnes. Most previous works, however, focus only on the dynamic behaviors of 
synchronization process in the stable region, i.e., global synchronization. In this letter, we demon- 
strate that synchronization process on complex networks can efficiently reveal the substructures of 
networks when the coupling strength of chaotic oscillators is under the lower boundary of stable 
region. Both analytic and numerical results show that the nodes belonging to the same component 
in the hierarchical network are tightly clustered according to the Euclidean distances between the 
state vectors of the corresponding oscillators, and different levels of hierarchy can be systematically 
unfolded by gradually tuning the coupling strength. When the coupling strengths exceed the upper 
boundary of stable region, the hierarchy of the network cannot be recognized by this approach. Ex- 
tensive simulations suggest that our method may provide a powerful tool to detect the hierarchical 
community structure of complex systems and deep insight into the relationship between structure 
and dynamics of complex systems. 
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Synchronization phenomena of interacting units, such 
as clapping peoples, fireflies and pendulums, have been 
noticed by scientists for a long time [l|-0|- For conve- 
nience of research, the interacting relations between units 
are usually abstracted as a network, whose nodes repre- 
sent units and edges indicate interactions between them. 
With the dramatic progress made in complex network 
science in the last decade, researchers from various dis- 
ciplines such as biology, physics and engineering com- 
munities have begun to explore a common interesting 
problem: the relation between topology of complex in- 
teractions and the emergent synchronization process 
[l3 |. On one hand, lots of works have investigated how 
synchronization process appears under different topology. 
For instance, the stability, one of the most important 
properties of synchronization, has been studied by the 
master stable function (MSF) method when the coupling 
networks have different topological properties like degree 
distribution, average path length, cluster coefficient and 
edge weight jl5l - [21| . In most cases, the stability is pro- 
portional to the eigenratio where Xmax is the 
largest eigenvalue of Laplacian matrix of network and A2 
is the second smallest one. Based on the studies of sta- 
bility, many methods are introduced to enhance the syn- 
chronizability via edge betweenness, topology modifica- 
tion, optimization, adaptive evolution, and so on [22M29j . 
On the other hand, few works focus on how topology of 
complex interactions can be refiected by synchronization 
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process. In [so], the hierarchical structure of network is 
gradually revealed from phase correlations in Kuramoto 
model [31] when the system evolves to the stable state. 
Recently, Ren et.al developed a universal approach to 
predict the exact network topology based solely on mea- 
suring the dynamical correlations of time series generated 
by the global synchronization [32|. 

In this letter, we focus on the problem of how sub- 
structures of complex networks can be refiected in the 
unstable synchronization processes taking place on net- 
works. Both the analytical and numerical results show 
that if the coupling strength is under the lower bound- 
ary of stable region, the nodes belonging to different 
components of hierarchical structure can be directly dis- 
tinguished by the state of the corresponding oscillator, 
which dramatically reduces the computational complex- 
ity of obtained cross-correlation among multiple time se- 
ries in traditional methods. Moreover, for a hierarchical 
network, we can unravel different levels of hierarchies in a 
top-down (or bottom-up) way by decreasing (or increas- 
ing) coupling strength gradually. However, if the cou- 
pling strength exceeds the upper boundary of the stable 
region, we cannot find the cluster structures of nodes even 
they couple in a hierarchical topology. In the following, 
we will first show that the states of coupling oscillators 
and substructure of network topology are both associ- 
ated with the Laplacian eigenvalues and eigenvectors in 
a theoretical way. 

Consider a generic system composed of N coupling 
chaotic oscillators, each one is represented by a m- 
dimensional state vector and ruled by the differential 
equations = F(x^). The network topology are defined 
by the Laplacian matrix L = {lij)NxN^ in which lij = 1 
if oscillator i and j are coupled, lij = otherwise, and 
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hi = ki (the degree of oscillator i). Therefore, we have 
the evolution equation of system 



N 



F(x,)-a^/,,H(x,). 



(1) 



Here we only consider undirected and binary networks, 

As the Laplacian matrix has zero row-sum, there is a 
global synchronization state of oscillators. 



Xi(t) X2(t) 



(2) 



Let (5x^ = x^ — s be the deviation of the ith oscillator 
from the synchronization manifold, we have 



N 



S±i = DF(s)(5x^ - crL)H(s) Gijd^j. 



(3) 



Let Ai < A2 < ... < Aat denote Laplacian eigenvalues 
and vi, V2, Vat denote the associated Laplacian eigen- 
vectors. To diagonalize the variational equations, we 
rewrite (5x as ^x = O^, where O = [vi, V2, vat] and 
i = ^2, •••5 ^at]"^- Hence, we obtain the N decoupled 
eigenmodes of variational equations 

ii = [Z)F(s) - (TA;iZ)H(s)]Ci, i = 1, 2, N, (4) 

Eq. dll is called the MSF of the complex system, 
is the coefficient corresponding to eigenvector in the 
linear combination. The eigenvector vi = [1,1,...,!] of 
Ai = corresponds to the synchronization manifold. The 
global synchronization of oscillators is achieved if aXi {i > 
2) totally falls in the stable region S which makes the 
maximum transverse Lyapunov exponents of Eq. ^ {i > 
2) be negative. The stable region S can be empty (^), 
bounded([ai a2]) or unbounded([ai oo]) with different 
chaotic oscillators and coupling function H. 

To identify different communities in a hierarchical net- 
work, we set coupling strength out of the stable region 
to ensure the differences between oscillator states remain. 
It's discovered that eigenvalues and eigenvectors of Lapla- 
cian matrix are strongly associated to substructures of 
network In a network with m communities, the 

former m — 1 non-trivial eigenvalues (from A2 to A^) are 
approximately equal to and much smaller than the oth- 
ers, which leads to a gap between the former m eigenval- 
ues and the others. And if a network is organized in a hi- 
erarchical way, there are more gaps between eigenvalues, 
which implies the different levels of hierarchy [30]. The 
eigenvectors of the former m — 1 eigenvalues describe the 
organization of network. Elements of eigenvectors corre- 
sponding to nodes in the same community share approx- 
imately the same value. In the hierarchical network, the 
eigenvectors corresponding to the eigenvalues before the 
ffist gap only show a coarse macroscopic organization of 
the highest hierarchy. When eigenvalues increasing and 
crossing following gaps, their eigenvectors can suggest a 
finer mesoscopic organization of lower hierarchies. We 
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FIG. 1: (a) eigenvalues Ai of Laplacian matrix of the hier- 
archical network considered. (b) Coupling strength am versus 
m. 



choose the coupling strength according to the inequality 
aXm < ai < crA^+i, which makes the coefficients ^2 to 
£^rn diverge so that the final states of oscillators are linear 
combinations of eigenvectors V2 to and the perturba- 
tions of v^+i to Vat vanish for sufficient evolving time. 
By increasing m, the more eigenvectors corresponding to 
larger eigenvalues are added as the basis of the linear 
combination of state vectors and more details of hierar- 
chical structure are uncovered by the state vectors. 

We simulate the synchronization process of coupling 
chaotic oscillators on a hierarchical network to conffim 
the above theoretical analysis. We ffistly choose the 
Rossler oscillators as an example of our investigation, 
whose stable region is bounded in [ai 0^2]. Their state 
vectors and differential equations can be written as x = 
[x, y, and F = [-{y + z), x + 0.2y, 0.2 + z{x - 9)f. 
Units are coupled by the linear coupling function H(x) = 
[x, 0, 0]-^. We therefore can get ai ~ 0.2, and a2 ^ 5 ac- 
cording to Eq. m The hierarchical network is ffistly 
proposed in [s^, which includes two hierarchical lev- 
els of components. Specifically, it consists 256 nodes, 
in which 16 non-overlapping clusters each containing 16 
nodes represent the first hierarchical level and 4 larger 
non-overlapping clusters each containing four ones from 
the ffist level consist the second hierarchical level. The 
connections of nodes at the ffist level (^ini), the sec- 
ond level {zin2) and with the rest nodes (zout) satisfy 
Zini + Zin2 + Zout = 18, and the hierarchical network is 
therefore indicated as Zini — ^in2- 

The synchronization process is performed on the 14 — 3 
hierarchical network. The coupling strength is defined as 

<Jra = ^[^ + T^],m = 2,3,...,N-l., (5) 
Z Am 

which guarantees aA^ < ai < aA^+i (i.e., the global 
synchronization is impossible). We show the Laplacian 
eigenvalues of the hierarchical network in Fig. [D^a), in 
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FIG. 2: (Color online) Distance matrix of 14-3 network with m changing from 2 to 19. With small m, only large clusters are 
recognized. Details of hierarchical structure emerge when the coupling strength decreases (i.e., m increases). 



which there are two obvious gaps that imply the two lev- 
els of hierarchy. As shown in Fig. [D^a), the first gap 
is between A4 and A5, while the second one is between 
A16 and A17. Thus, we are able to unfold the hierarchi- 
cal structure gradually as m increases from 2 to 19. By 
considering the eigenvalues and lower limit of the stable 
region of Rossler oscillator, we set the coupling strength 
according to Eq. O which monotonously decreases with 
m (see in Fig. [Hb)). Moreover, the maximum coupling 
strength (72 is also constrained by (72 A at < 0^2 to make the 
coefficient corresponding to Xi{i > m + 1) converge and 
perturbations of v^+i to vat vanish for sufficient evolv- 
ing time. 

We use the Euler method with time step At = 0.001 
and total evolving time T = 100 to numerically simulate 
the synchronization process. The initial states of Rossler 
oscillators (i.e., the nodes of hierarchical network) are 
uniformly distributed in the interval [0, 1]. The difference 
of state between the nodes is suggested by their Euclidean 
distance dij 



{Xi - Xj)^ + {yi - VjY + {Zi 

Figure. [2] shows the distance matrices at T = 100 at 
a changing coupling strength cr^ from m = 2tom = 19. 
Each distance matrix is averaged over a hundred runs of 
synchronization process. We find that the blocks consist- 
ing of trivial values exist in all distance matrices. Inside 
each block, the distances between nodes are small, while 
the distance between nodes across different blocks is large 
enough to unfold the hierarchy of network. Specifically, 
with smaller m, the coupling strength is large, which 
causes only the coefficients corresponding to the small- 
est eigenvalues to diverge. In these cases, we can only rec- 
ognize the larger clusters at the second hierarchical level. 



When m > 5, the coefficients corresponding to the eigen- 
values after the ffist gap begin to diverge, which suggests 
that in the distance matrices, the 4 non-overlapping com- 
ponents at the second hierarchical level split into smaller 
components at the ffist hierarchical level. As m increases 
further from 5 to 16, the eigenvectors corresponding to 
the larger eigenvalues can provide the difference between 
nodes that belong to the same component at the second 
hierarchical level, as is manifested in Fig [2l Thus, the 
16 non-overlapping components at the second hierarchi- 
cal level become much more obvious and the two-level 
hierarchy of network is gradually revealed. For m > 5, 
the coupling strength decrease to trivial values so that 
the synchronization process is weakened, which can help 
reveal more detailed structure of hierarchical network. 

To show the clustering of nodes more explicitly, we 
present the dendrogram plots of agglomerative hierarchi- 
cal trees and hierarchical structures implied by the hi- 
erarchical trees (see Fig. [3]). The dendrogram plots of 
agglomerative hierarchical trees for m = 2,3,4 and 16 
are shown in Fig. [3] (a), (b), (c) and (d), respectively, 
and the corresponding hierarchical structures of network 
are described in Fig. [3](e), (f), (g) and (h). When m = 2, 
the eigenvector V2 is dominant and the distances between 
nodes is to a large extent determined by it. The network 
in this case is partitioned into two asymmetric clusters, 
and one is much larger than the other one (corresponding 
to the distance matrix of m = 2 in Fig. [2j). Furthermore, 
when m increases from 2 to 4, the four non-overlapping 
components at the second level of hierarchical network 
emerge one by one. At last, the 16 non-overlapping com- 
ponents at the ffist level and 4 non-overlapping compo- 
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FIG. 3: (Color online) Plots ((a), (b), (c) and (d)) show the 
hierarchical tree, and plots ((e), (f), (g) and (h)) describe the 
cluster structure of networks implied by the hierarchical tree 
for m = 2, 3, 4 and 16 respectively. It's demonstrated clearly 
that the hierarchical structure unfolded by synchronizing pro- 
cess change as coupling strength decreasing. 



nents at second level are totally unfolded as the coeffi- 
cients corresponding to eigenvalues from V2 to all 
diverge, as shown in Fig. ^d) and (h). 

We have demonstrated how the synchronization pro- 
cess of coupling chaotic oscillators can unravel the hier- 
archical structure of network at different levels by tuning 
the coupling strengths under the lower boundary of sta- 
ble region. On the other hand, the global synchroniza- 
tion of Rossler oscillators is also broken if the coupling 
strengths exceed the upper boundary of stable region. 
In this situation, states vectors are linear combinations 
of eigenvectors corresponding to the largest eigenvalues, 
of which the elements are almost randomly distributed 
rather than clustered according to the substructures of 
network. Thus, the hierarchical structure of network can- 
not be revealed by the distance matrices. In the simula- 
tion, we set the coupling strength in agreement with the 
the inequality crA^ < 0^2 < crA^+i(m = 2, 3, ...A^— 1). To 
eliminate the influence from the eigenvectors correspond- 
ing to small eigenvalues, the coupling strength is set to 
be larger than ai/X2 — 0.17. For instance, if the cou- 
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FIG. 4: (Color online) Distance matrix with coupling 
strength a = 0.19 at T = 100. No hierarchical structure 
can be seen in this matrix. 
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FIG. 5: (Color online) Distance matrix with coupling 
strength cr = 4, 3, 2 and 1 at T = 100. As coupling strength 
decreasing, different levels of hierarchical structures are grad- 
ually revealed by synchronizing process 



pling strength is chosen as a = 0.19, the smallest value 
of m that satisfies aXm < 0^2 < crAm+i is m = 229, which 
suggests that with i > 229 diverges. The synchroniza- 
tion process performed on the same hierarchical network 
also evolves for time T = 100, and the distance matrix 
is also averaged over a hundred realizations. As shown 
in Fig. m the visualization of distance matrix shows no 
evidence of clustering of nodes. Most values in distances 
matrix are trivial, because almost all the states vectors 
are synchronized to the same state. Only few nodes are 
dramatically far away from the others. Therefore, al- 
though the global synchronization is also impossible be- 
cause coupling strength exceeds the upper boundary of 
stable region, it does not provide useful information to 
unfold the hierarchical structure of networks. 

We carried out the comprehensive numerical experi- 
ments using Lorenz oscillators to demonstrate the gener- 
ality of our theoretical analysis. The state vector and dif- 
ferential equations of Lorenz oscillators can be written as 



5 



X = [x, y, z]^ and F = [I0{y — x), x(28 — z) — y,xy — 
respectively. The units are also coupled by the linear 
coupling function H(x) = [x,0,0]^. The stable region 
of linear coupling Lorenz oscillators is [aioo]. Note that 
the numeric method gives ai = 2, but the accurate value 
is much lareger in the numeric simulation of synchro- 
nization process, and we therefore don't set the coupling 
strength according to Eq. [5l Furthermore, because the 
upper boundary of stable region of Lorenz oscillators is 
infinite (oo), we only need to consider the low boundary 
of the stable region, and vary the coupling strength a 
from 4 to 1 for instance. The synchronization processes 
are performed on the 14-3 hierarchical network with the 
same conditions used in coupling Rossler oscillators. The 
distance matrices of nodes at T = 100 are shown in Fig. 
[5l We can see that the coupling Lorenz oscillators on the 
hierarchical network beyond global synchronization re- 
gion can gradually unfold the substructures and different 
hierarchical levels of the network. 

In conclusion, we have both theoretically and numeri- 
cally investigated the synchronization processes on hier- 
archical networks beyond stable region (i.e., the coupling 
strength is under the low boundary ai/A2), which can 
be used to unfold the fine substructure such as differ- 
ent levels of hierarchies in an efficient way. The results 



show that the nodes belonging to the same component 
are tightly clustered according to the state-vector dis- 
tances between the corresponding coupling oscillators. 
We also find that if the coupling strength exceeds the 
upper boundary 0^2 /Aat, the hierarchical structure of net- 
work cannot be effectively unfolded as the trajectories of 
the oscillators mix together. Since the nodes belonging to 
the same community can be directly identified by their 
state vectors, our approach can detect the hierarchical 
community structure of complex systems with very low 
computational complexity, compared to traditional cross- 
correlation based or spectral approaches. Our approach 
to unravel the hierarchical structure of network suggests 
that the dynamics of the interaction contains abundant 
information about the topology of the interaction, which 
can be utilized to infer the fine structure of complex sys- 
tem. 

This work is supported by the National Natural Sci- 
ence Foundation of China under Grant Nos. 60874090, 
60974079, and 61004102. S-MC appreciates the financial 
support of the Fundamental Research Funds for the Cen- 
tral Universities (No.WK2100230004). JZ acknowledges 
the financial support from the National Natural Natural 
Science Foundation of China under grant Nos. 61004104 
and 61104143. 



[1] C. Huygens, Horoloquium Oscilatorium (Apud F. 

Muguet, Paris, 1673). 
[2] L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. 64, 821 

(1990) 

[3] Z Neda, E Ravasz, Y Brechet, T. Vicsek, and A. L. 
Barabasi, Nature 403, 849 (2000). 

[4] A. Pikovsky, M. Rosenblum, and J. Kurths, Syn- 
chronization: A Universal Concept in Nonlinear Sci- 
ence (Cambridge University Press, New York, 2002). 

[5] S. H. Strogatz, Nature 410, 268 (2001). 

[6] E. Oh, K. Rho, and B. Kahng, Phys. Rev. E 72, 047101 
(2005). 

[7] C. S. Zhou, L. Zemanova, G. Zamora, C. Hilgetag, and 

J. Kurth, Phys. Rev. Lett. 97, 238103 (2006). 
[8] J. Gomez- Gar dries, Y. Moreno, and A. Arenas, Phys. 

Rev. Lett. 98, 034101 (2007) 
[9] Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno, and 

C.Zhou, Phys. Rep. 469, 93 (2008). 
[10] W. YU, G. Chen, and J. Lii, Automatica 45, 429 (2009). 
[11] T. Nishikawa and A. E. Motter, Proc. Nat. Acad. Sci. 

107, 10342 (2010). 
[12] Z. Zhuo, S. M. Cai, Z. Q. Fu, and W. X. Wang, New J. 

Phys. 13, 053030 (2011). 
[13] J. Zhang, K. Zhang, X. Xu, C. K. Tse and M. Small, 

New J. Phys. 11, 113003 (2009). 
[14] J. Zhang, C. Zhou, X. Xu, and M. Small, Phys. Rev. E 

82, 026116 (2010). 
[15] L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. 80, 2109 

(1998). 

[16] M. Barahona, L.M. Pecora, Synchronization in small- 
world systems, Phys. Rev. Lett. 89, 054101 (2002). 
[17] T. Nishikawa, A. E. Motter, Y. C. Lai, and F. C. Hop- 



pensteadt, Phys. Rev. Lett. 91, 014101 (2003). 
[18] H. Hong, B. J. Kim, M. Y. Choi, and H. Park, Phys. Rev. 

E 69, 067105 (2004). 
[19] C. Li and G. Chen, Physica A 341, 73 (2004). 
[20] S. H. Yook and H. Meyer- Ortmanns, Physica A 371, 781 

(2006) 

[21] L V. Belykh, V. N. Belykh, and M. Hasler, Physica D 

195, 188 (2004). 
[22] M. Chavez, D. U. Hwang, A. Amann, H. G. E. Hentschel, 

and S. Boccaletti, Phys. Rev. Lett. 94, 218701 (2005). 
[23] M. Zhao, T. Zhou, B. H. Wang, and W. X. Wang, Phys. 

Rev. E 72, 057102 (2005). 
[24] C. Y. Yin, W. X. Wang, G. Chen, and B. H. Wang, Phys. 

Rev. E 74, 047102 (2006). 
[25] L. Donetti, P. L Hurtado, and M. A. Mufioz, Phys. Rev. 

Lett. 95, 188701 (2005). 
[26] T. Nishikawa and A. E. Motter, Physica D 224, 77 

(2006). 

[27] T. Nishikawa and A. E. Motter, Phys. Rev. E 73, 

065106(R) (2006). 
[28] C. Zhou and J. Kurths, Phys. Rev. Lett. 96, 164102 

(2006). 

[29] G. Yan, G. Chen, J. H. Lii, and Z. Q. Fu, Phys. Rev. E 
80, 056116 (2009). 

[30] A. Arenas, A. Diaz-Guilera, and C. J. Perez- Vicente, 
Phys. Rev. Lett. 96, 114102 (2006). 

[31] Y. Kuramoto, Chemical Oscillations, Waves, and Turbu- 
lence (Springer, Berlin, 1984). 

[32] J. Ren, W. X. Wang, B. W. Li, and Y. C. Lai, Phys. Rev. 
Lett. 104, 058701 (2010) 

[33] A. Capocci, V. D. P. Servedio, G. Galdarelli, and F. Co- 
laiori, Physica A 352, 669 (2005). 



6 



[34] L. Huang, K. Park, Y. C. Lai, L. Yang, and K. Q. Yang, 
Phys. Rev. Lett. 97, 164101 (2006) 



